home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / raytrace / radiance / simplerd.lha / simplerad / FinalFTP / Light / ff.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-05-21  |  26.7 KB  |  792 lines
  1. /**********************************************************************/
  2. /* ff.c                                                               */
  3. /*                                                                    */
  4. /* Routines for form factor approximation. Algorithms used:           */
  5. /*                                                                    */
  6. /* Wallace89 : differential area to disc approximation, with uniform  */
  7. /*             or jittered sampling, and adaptive source sampling.    */
  8. /*             Change from circle to ellipse approximation, plus add  */
  9. /*             in differential to finite approximation for            */
  10. /*             perpendicular case                                     */
  11. /* Baum/Rush90 : analytic ff approximation for special cases:         */
  12. /*              1) Proximity violations.                              */
  13. /*              2) Visibility violations.                             */
  14. /*              3) PR violations.                                     */
  15. /*              4) Wallace89 problems.                                */
  16. /* Hanrahan91 : BF refinement, and O(n) form-factor interaction       */
  17. /* (not completed in this revision).                                  */
  18. /*                                                                    */
  19. /* Copyright (C) 1992, Bernard Kwok                                   */
  20. /* All rights reserved.                                               */
  21. /* Revision 1.0                                                       */
  22. /* May, 1992                                                          */
  23. /**********************************************************************/
  24. #include <stdio.h>
  25. #include <stdlib.h>
  26. #include <math.h>
  27. #include "geo.h"
  28. #include "misc.h"
  29. #include "io.h"
  30. #include "struct.h"
  31. #include "rad.h"
  32. #include "ray.h"
  33. #include "ff.h"
  34. #include "rtime.h"
  35.  
  36. /**********************************************************************/
  37. extern OptionType Option;
  38. extern char *ProgName;
  39. extern int HitObject();
  40. extern Time_Stats tstats;
  41.  
  42. /**********************************************************************/
  43. /* Form factor globals                                                */
  44. /**********************************************************************/
  45. double MAX_F_DIFF_EDGE = 0.016;       /* maximum allowed difference 
  46.                      between ff of 2 vertices of a 
  47.                      edge before you have to 
  48.                      subdivide */
  49. double MIN_ELEMENT_PERC = 0.1;        /* minimum % area of element    
  50.                      for subdivision */
  51.  
  52. FF_OptionType FF_Options = {          /* Default B, F tolerances and */
  53.   DEFAULT_B_min, DEFAULT_B_max,       /* default uniform sampling    */
  54.   DEFAULT_F_min, DEFAULT_F_max,       /* with default 4x4 sampling   */
  55.   UNIFORM, DEFAULT_SAMPLES,           /* grid */
  56.   NUMERICAL_FF, /* MIN_ELEMENT_AREA, */
  57.   MAX_F_DIFF_POLY, /* MAX_F_DIFF_EDGE */ 0.016,
  58.   MAX_PATCH_LEVELS, DISC_FF
  59.   };
  60.  
  61. /**********************************************************************/
  62. /* Print form factor options                                          */
  63. /**********************************************************************/
  64. void ff_print_FF_Options(ffopt)
  65.      FF_OptionType ffopt;
  66. {
  67.   FILE *fp;
  68.   if (Option.device == PRINT)
  69.     fp = stdout;
  70.   else
  71.     fp = Option.StatFile;
  72.  
  73.   fprintf(fp,"\n\t*** Form-factor options *** \n");
  74.   fprintf(fp,"\tB {max %g} {min %g}\n", ffopt.B_max, ffopt.B_min);
  75.   fprintf(fp,"\tF {max %g} {min %g}\n", ffopt.F_max, ffopt.F_min);
  76.   fprintf(fp,"\t%s %s form factors", 
  77.      (ffopt.fftype == ANALYTIC_FF ? "Analytic" : "Numerical"),
  78.      (ffopt.sample_shape == ELLIPSE_FF ? "elliptical" : "circular"));
  79.   if (ffopt.use_analytic) fprintf(fp,", WITH Analytic for special cases.\n");
  80.   else fprintf(fp,".\n");
  81.   fprintf(fp,"\tMinimum subdivision area %g\n", ffopt.min_element_area);
  82.   fprintf(fp,"\tMaximum poly form factor difference %g\n", ffopt.F_diff_poly);
  83.   fprintf(fp,"\tMaximum edge form factor difference %g\n", ffopt.F_diff_edge);
  84.   fprintf(fp,"\tMaximum patch subdivision levels %d\n", ffopt.max_levels);
  85.  
  86.   if (Option.ff_raytrace) {
  87.     fprintf(fp,"\tUsing ray cast visibility testing, with\n");
  88.     fprintf(fp,"\t%s sampling with %d %s samples per polygon.\n", 
  89.        (ffopt.sampling_type == UNIFORM ? "uniform" : "jittered"),
  90.        ffopt.num_samples,
  91.        (ffopt.varying_numsamps == 1 ? "adjusted" : "unadjusted"));
  92.   }
  93. }
  94.   
  95. /**********************************************************************/
  96. /* Compute xyz coordinates of a point on quadrilateral given it's u,v */
  97. /* coordinates using bi-linear mapping                                */
  98. /**********************************************************************/
  99. void UVToXYZ(quad, u, v, xyz)
  100.      Vector quad[MAX_PATCH_VTX];
  101.      double u,v;
  102.      Vector *xyz;
  103. {
  104.   xyz->x = quad[0].x * (1.-u)*(1.-v) + quad[1].x * (1.-u)*v + 
  105.     quad[2].x *u*v + quad[3].x * u*(1.-v);
  106.   xyz->y = quad[0].y * (1.-u)*(1.-v) + quad[1].y * (1.-u)*v + 
  107.     quad[2].y *u*v + quad[3].y * u*(1.-v);
  108.   xyz->z = quad[0].z * (1.-u)*(1.-v) + quad[1].z * (1.-u)*v + 
  109.     quad[2].z *u*v + quad[3].z * u*(1.-v);
  110. }
  111.  
  112. /**********************************************************************/
  113. /* Initialize the set of samples */
  114. /**********************************************************************/
  115. int Init_Samples(s)
  116.      Sampletype s[MAX_SAMPLES];
  117. {
  118.   int i;
  119.   for(i=0;i<MAX_SAMPLES;i++) { s[i].p = 0; }
  120. }
  121.  
  122. /**********************************************************************/
  123. /* Print set of samples (for testing only)                            */
  124. /**********************************************************************/
  125. void ff_print_samples(s,num_samples,label)
  126.      Sampletype s[MAX_SAMPLES];
  127.      int num_samples;
  128.      char *label;
  129. {
  130.   int i,j;
  131.  
  132.   printf("\n\tSamples %s (%d) {\n", label, num_samples);
  133.   for(i=0;i<num_samples;i++) {
  134.     printf("\tPoly[%d] %s%d {", i, s[i].p->name, s[i].p->id);
  135.     printf(" Pos %g %g %g;", s[i].pos.x, s[i].pos.y, s[i].pos.z);
  136.     printf(" Norm %g %g %g;", s[i].norm.x, s[i].norm.y, s[i].norm.z);
  137.     printf(" A %g\n", s[i].area);
  138.  
  139.     printf("\tCorners: ");
  140.     for(j=0;j<s[i].p->numVert;j++)
  141.       printf("(%g %g %g) ", s[i].corners[j].x, s[i].corners[j].y, 
  142.          s[i].corners[j].z);
  143.     printf("\n");  
  144.   }
  145.   printf("\t}\n");  
  146. }
  147.  
  148. /**********************************************************************/
  149. /* Sample a triangle's centre, and recursively subdivide uniformly */
  150. /**********************************************************************/
  151. void Sample_Triangle_Center(tri,s,iter,sample_num)
  152.      Vector tri[4];
  153.      Sampletype s[MAX_SAMPLES];
  154.      int iter;
  155.      int *sample_num;
  156. {
  157.   int i;
  158.   Vector mid[3];
  159.   Vector corners[3];
  160.   Vector tri_sample;
  161.   
  162.   if (iter >= 0) 
  163.     iter--;
  164.   else
  165.     return;
  166.  
  167.   /* Sample center of triangle, store corners of sample */
  168.   if (iter == -1) {
  169.     tri_sample.x = (tri[0].x + tri[1].x + tri[2].x) / 3.0;
  170.     tri_sample.y = (tri[0].y + tri[1].y + tri[2].y) / 3.0;
  171.     tri_sample.z = (tri[0].z + tri[1].z + tri[2].z) / 3.0;
  172.     s[*sample_num].pos = tri_sample;
  173.     for(i=0;i<3;i++) 
  174.       s[*sample_num].corners[i] = tri[i];
  175.     (*sample_num)++;
  176.   }
  177.  
  178.   /* Find middle of edges of triangle */
  179.   mid[0] = *vmiddle(&tri[0],&tri[1]);
  180.   mid[1] = *vmiddle(&tri[1],&tri[2]);
  181.   mid[2] = *vmiddle(&tri[2],&tri[0]);
  182.  
  183.   /* Recursively sample - check only triangles not checked before */
  184.   corners[0] = tri[0];  corners[1] = mid[0];  corners[2] = mid[2];
  185.   Sample_Triangle_Center(corners,s,iter,sample_num);  
  186.  
  187.   corners[0] = mid[0];  corners[1] = tri[1];  corners[2] = mid[1];
  188.   Sample_Triangle_Center(corners,s,iter,sample_num);
  189.  
  190.   corners[0] = mid[2];  corners[1] = mid[0];  corners[2] = mid[1];
  191.   Sample_Triangle_Center(corners,s,iter,sample_num);  
  192.  
  193.   corners[0] = mid[2];  corners[1] = mid[1];  corners[2] = tri[2];
  194.   Sample_Triangle_Center(corners,s,iter,sample_num);  
  195.  
  196. }
  197.  
  198. /**********************************************************************/
  199. /* Sample subdivided triangle */
  200. /**********************************************************************/
  201. void Sample_Triangle(tri,s,num_samples,sample_type, corners, area, normal)
  202.      Polygon *tri;
  203.      Sampletype s[MAX_SAMPLES];
  204.      int num_samples;
  205.      int sample_type;
  206.      Vector corners[MAX_PATCH_VTX];
  207.      double area;
  208.      Vector normal;
  209. {
  210.   int i;
  211.   int sample_num = 0;
  212.   int iter;
  213.   double result;
  214.  
  215.   /* Store rest of info first */
  216.   for(i=0;i<num_samples;i++) {
  217.     s[sample_num].p = tri;
  218.     s[sample_num].norm = normal;
  219.     s[sample_num++].area = area / num_samples;
  220.   }
  221.  
  222.   /* Find sample positions, by recursively subdividing the triangle for
  223.      "iter" number of iterations */
  224.   iter = 0; result = num_samples / 4.0;
  225.   while (result >= 1.0) {
  226.     iter++;
  227.     result /= 4.0;
  228.   }
  229.   sample_num = 0;
  230.   Sample_Triangle_Center(corners,s,iter,&sample_num);
  231. }
  232.   
  233. /**********************************************************************/
  234. /* Sample quadralateral */
  235. /**********************************************************************/
  236. void Sample_Quad(quad,s,num_samples,sample_type, corners, area, normal)
  237.      Polygon *quad;
  238.      Sampletype s[MAX_SAMPLES];
  239.      int num_samples;
  240.      int sample_type;
  241.      Vector corners[MAX_PATCH_VTX];
  242.      double area;
  243.      Vector normal;
  244. {
  245.   int i,j;
  246.   double u,v, du, dv;
  247.   double jitter_u, jitter_v;
  248.   Vector tmp;
  249.   int nu, nv;
  250.   int sample_num = 0;
  251.  
  252.   for(i=0;i<MAX_PATCH_VTX;i++) 
  253.     corners[i] = quad->vert[i]->pos;
  254.   
  255.   /* Set grid subdivision */
  256.   nv = nu = CEILING(sqrt((double)num_samples));
  257.   dv = du = 1.0 / (double) nu;
  258.   
  259.   for(i=0,u=du/2.0; i<nu; i++, u+=du) {
  260.     for(j=0,v=dv/2.0; j<nv; j++, v+=dv, sample_num++) {
  261.       UVToXYZ(corners, u, v, &tmp);
  262.       s[sample_num].pos = tmp;
  263.       UVToXYZ(corners, u-du/2.0, v-dv/2.0, &tmp);
  264.       s[sample_num].corners[0] = tmp;
  265.       UVToXYZ(corners, u+du/2.0, v-dv/2.0, &tmp);
  266.       s[sample_num].corners[1] = tmp;
  267.       UVToXYZ(corners, u+du/2.0, v+dv/2.0, &tmp);
  268.       s[sample_num].corners[2] = tmp;
  269.       UVToXYZ(corners, u-du/2.0, v+dv/2.0, &tmp);
  270.       s[sample_num].corners[3] = tmp;
  271.       s[sample_num].p = quad;
  272.       s[sample_num].norm = normal;
  273.       s[sample_num].area = area / num_samples;
  274.     }
  275.   }
  276. }
  277.  
  278. #define MIN_SOURCE_UNSHOT 0.1
  279. #define MAX_SOURCE_UNSHOT 0.5
  280. #define USE_NEW 1
  281. /**********************************************************************/
  282. /* Routine to adjust number of samples to generate based on given     */
  283. /* criteria. Currently adjust # of rays depending on:                 */
  284. /*   - level of subdivsion of element.                                */
  285. /*   - closeness and area of light wrt receiver (i.e. what solid      */
  286. /*     angle does it subtend?)                                        */
  287. /*   - total radiosity of source < tolerance                          */
  288. /**********************************************************************/
  289. int Adjusted_NumSamples(p)
  290.      Polygon *p;
  291. {   
  292.   int i;
  293.   int num_samps, scale;
  294.   int sample_weight; /* 0 = 1 sample, 1 = 4 samples, 2 = sixteen samples */
  295.  
  296. #ifdef USE_NEW
  297.   if (FF_Options.num_samples == 16)
  298.     sample_weight = 2;
  299.   else if (FF_Options.num_samples == 4)
  300.     sample_weight = 1;
  301.   else 
  302.     sample_weight = 0;
  303.  
  304.   /* Adjust based on source radiosity */
  305.   for (i=0;i<MAX_SPECTRA_SAMPLES;i++)
  306.     if (p->unshot_B.samples[i] < MIN_SOURCE_UNSHOT) 
  307.       sample_weight--;
  308.  
  309.   /* Adjust based on solid angle source subtends receiver */
  310.   /* (Not implemented) */
  311.  
  312.   /* Adjust based on receiver size -- new */
  313.   if (p->level > 0) {
  314.     scale = (int) Power(4.0,(double)p->level);
  315.     if (scale >= FF_Options.num_samples)
  316.       sample_weight = 0;
  317.     else
  318.       sample_weight -= p->level;
  319.   }
  320.  
  321.   num_samps = 1;
  322.   if (sample_weight < 0)      sample_weight = 0;
  323.   else if (sample_weight > 2) sample_weight = 2;
  324.   for (i=1;i<sample_weight;i++)
  325.     num_samps = num_samps * 4;
  326.  
  327.   if (Option.statistics) {
  328.     if (Option.device == PRINT)
  329.       printf("\tAdjusted number samples = %d\n", num_samps);
  330.     else
  331.       fprintf(Option.StatFile,"\tAdjusted number samples = %d\n", num_samps);
  332.   }
  333.  
  334. #else
  335.   /* Adjust based on receiver size -- old */
  336.   if (p->level > 0) {
  337.     scale = (int) Power(4.0,(double)p->level);
  338.     if (scale >= FF_Options.num_samples)
  339.       num_samps = 1;
  340.     else
  341.       num_samps = FF_Options.num_samples /  scale;
  342.   } else 
  343.     num_samps = FF_Options.num_samples;
  344.  
  345. #endif
  346.  
  347.   return num_samps;
  348. }
  349.  
  350. /**********************************************************************/
  351. /* Generate samples (s) on a patch (p)                                */
  352. /**********************************************************************/
  353. int ff_Generate_Patch_Samples(p,s,corners,area,normal)
  354.      Polygon *p;
  355.      Sampletype s[MAX_SAMPLES];
  356.      Vector corners[MAX_PATCH_VTX];
  357.      double area;
  358.      Vector normal;
  359. {
  360.   int num_samples;
  361.   int i;
  362.  
  363.   if (FF_Options.fftype == ANALYTIC_FF)
  364.     num_samples = DEFAULT_ASAMPLES;
  365.   else {
  366.     if (FF_Options.varying_numsamps)
  367.       num_samples = Adjusted_NumSamples(p);
  368.     else
  369.       num_samples = FF_Options.num_samples;
  370.   }
  371.  
  372.   switch (FF_Options.sampling_type) {
  373.  
  374.   case UNITEST:               /* Test only */
  375.  
  376.     for (i=0;i<num_samples;i++) {
  377.       s[i].p = p;
  378.       s[i].pos = corners[i];
  379.       s[i].norm = normal;
  380.       s[i].area = area / num_samples;
  381.     }
  382.     break;
  383.     
  384.   case UNIFORM:               /* Sample centres of grid patch */
  385.  
  386.     if (p->class == PATCH) 
  387.       Sample_Quad(p,s,num_samples,UNIFORM,corners,area,normal);
  388.     else if (p->class == TRIANGLE)
  389.       Sample_Triangle(p,s,num_samples,UNIFORM,corners,area,normal);
  390.     break;
  391.     
  392.   case JITTERED:                /* Sample with centers jittered */
  393.  
  394.     if (p->class == PATCH) 
  395.       Sample_Quad(p,s,num_samples,JITTERED,corners,area,normal);
  396.     else if (p->class == TRIANGLE)
  397.       Sample_Triangle(p,s,num_samples,JITTERED,corners,area,normal);
  398.     break;
  399.     
  400.   default:
  401.     fprintf(stderr,"%s: Unknown sampling type %d\n", ProgName,
  402.         FF_Options.sampling_type);
  403.     exit(1);
  404.   }
  405.   return (num_samples);
  406. }
  407.  
  408.  
  409. /**********************************************************************/
  410. /* Test if a reciever vertex is on the plane of the src poly */
  411. /**********************************************************************/
  412. int Patches_Abut(src,rec)
  413.      Polygon *src, *rec;
  414. {
  415.   int i;
  416.   Vector rv, sv, sn, ray_dir;
  417.   double cos_angle_p;
  418.   int is_perp = FALSE;
  419.   int is_touching = FALSE;
  420.  
  421.   /* Is perpendicular test */
  422.   sn = src->normal[0];
  423.   i = 0;
  424.   while (i<rec->numVert && is_perp == FALSE) {
  425.     sv = src->vert[0]->pos;
  426.     rv = rec->vert[i]->pos;
  427.     ray_dir = *vsub(&sv,&rv);
  428.     norm(&ray_dir);
  429.     cos_angle_p = dot(&sn, vnegate(&ray_dir));
  430.     if (cos_angle_p == 0.0) 
  431.       is_perp = TRUE;
  432.     i++;
  433.   }
  434.  
  435.   if (is_perp) {
  436.     i = 0;
  437.     while (i<src->numVert && is_touching == FALSE) {
  438.       sv = src->vert[i]->pos;
  439.       if (vsame(&rv, &sv, DAMN_SMALL)) is_touching = TRUE;
  440.       i++;
  441.     }
  442.   }
  443.   return is_touching;
  444. }  
  445.  
  446. int FF_Testing = 0;                               /* Just for testing */
  447. /**********************************************************************/
  448. /* Shoot a ray from vertex to sample point on source                  */
  449. /* Return intersection distance                                       */
  450. /**********************************************************************/
  451. int shoot_ray(ray_orig, ray_dir, s, t)
  452.      Vector ray_orig, ray_dir;
  453.      Sampletype s;
  454.      double *t;
  455. {
  456.   double tmp;
  457.   Vector normray_dir;
  458.  
  459.   tmp = vlen(&ray_dir);      /* Find distance from source to receiver */
  460.   *t = tmp;
  461.   if (FF_Testing) {
  462.     printf("\tShooting from %g,%g,%g to %g,%g,%g ", 
  463.        ray_orig.x, ray_orig.y, ray_orig.z,
  464.        ray_dir.x, ray_dir.y, ray_dir.z);
  465.     printf("at pos %g %g %g. t = %g\n", s.pos.x, s.pos.y, s.pos.z, *t);
  466.   }
  467.  
  468.   /* Find out if ray hit any object before reaching sample point on
  469.      source */
  470.   if (Option.visibility == FORM_FACTOR) {
  471.     normray_dir = ray_dir;
  472.     norm(&normray_dir);
  473.     if (HitObject(ray_orig, normray_dir, *t))
  474.       return 0;
  475.     else
  476.       return 1;
  477.   } else
  478.     return 1;
  479. }
  480.  
  481. /*====================================================================*/
  482. /* Analytically integrated form-factor from [BAUM90], using contour   */
  483. /* integration.                                                       */
  484. /*====================================================================*/
  485. double ff_da_analytic(pptr, from, corners, pnorm)
  486.      Polygon *pptr;                    /* Source polygon */
  487.      Vector *from;                     /* Receiver vertex */
  488.      Vector *pnorm;                    /* Receiver normal */
  489.      Vector corners[MAX_PATCH_VTX];    /* Source sample corners */
  490. {
  491.   int i, ii;
  492.   double dff;
  493.   Vector R[MAX_PATCH_VTX];
  494.   Vector Tau[MAX_PATCH_VTX];
  495.   Vector direction;
  496.   double corner_angle;
  497.   double cos_corner_angle;
  498.  
  499.   /* Set receiver normal, and rays from receiver 
  500.      position to corners of source */
  501.   for(i=0;i<(pptr->numVert);i++) {
  502.     R[i] = *vsub(&corners[i], from);
  503.     norm(&R[i]);
  504.   }
  505.   
  506.   /* Calculate ff by numerically integrating */
  507.   dff = 0.0;
  508.   for(i=0;i<pptr->numVert;i++) {
  509.     ii = (i + 1) % pptr->numVert;
  510.     cos_corner_angle = dot(&R[i], &R[ii]);
  511.     corner_angle = acos(cos_corner_angle);
  512.     direction = cross(&R[i], &R[ii]);
  513.     norm(&direction);
  514.     Tau[i] = *vscalar(corner_angle, &direction);
  515.     dff += fabs(dot(pnorm, &Tau[i]));
  516.   }
  517.  
  518.   dff = dff / (2.0 * PI);
  519.   return(dff);
  520. }
  521.  
  522. #define MAX_SRC_TO_RECVTX_ENERGY 5.0 
  523. /*====================================================================*/
  524. /* Ray tracing ff -- circular disc approximation by [Wallace89]       */
  525. /*   1) Assume uniform weighted area per sample                       */
  526. /*   2) Uniform or jittered sample points on area samples             */
  527. /*====================================================================*/
  528. double ff_da_disc(v,n,p,samples,num_samples)
  529.      Vector v, n;                /* Vertex and normal to sample from */
  530.      Polygon *p;                       /* Polygon to be sampled      */
  531.      Sampletype samples[MAX_SAMPLES];  /* Samples on source          */
  532.      int num_samples;                  /* Number of source samples   */
  533. {
  534.   int i;
  535.   double ff, cos_angle_v, cos_angle_p, t, dff, w; /* for ff calc     */
  536.   Vector ray_dir;                    /* Direction of ray from vtx    */
  537.   Vector norm_ray_dir;               /* Normalized ray direction     */
  538.                     
  539.   /* Calculate delta form-factor from each sample */
  540.   ff = 0.0;
  541.   for (i=0;i<num_samples;i++) {
  542.     ray_dir = *vsub(&(samples[i].pos), &v);
  543.     norm_ray_dir = ray_dir;
  544.     norm(&norm_ray_dir);
  545.     w = samples[i].area;
  546.     if (shoot_ray(v, ray_dir, samples[i], &t)) {
  547.       cos_angle_v = dot(&n, &norm_ray_dir);
  548.       cos_angle_p = dot(&(samples[i].norm), vnegate(&norm_ray_dir));
  549.       if (cos_angle_p == 0.0) cos_angle_p = cos(DTOR * 90.0);      
  550.       dff = (cos_angle_v * cos_angle_p) / ((M_PI * t*t) + samples[i].area);
  551.       
  552.       /* Take area weighted average of samples as form factor */
  553.       dff = dff * w;
  554.       ff += dff;
  555.     }
  556.   }
  557.   return(ff);
  558. }
  559.  
  560. /*====================================================================*/
  561. /* Ray tracing ff -- elliptical disc approximation [Siegal81]         */
  562. /*====================================================================*/
  563. double ff_da_ellipse(v,n,p,samples,num_samples)
  564.      Vector v, n;                /* Vertex and normal to sample from */
  565.      Polygon *p;                       /* Polygon to be sampled      */
  566.      Sampletype samples[MAX_SAMPLES];  /* Samples on source          */
  567.      int num_samples;                  /* Number of source samples   */
  568. {
  569.   int i;
  570.   double ff, cos_angle_v, cos_angle_p, t, dff, w; /* for ff calc     */
  571.   Vector ray_dir;                    /* Direction of ray from vtx    */
  572.   Vector norm_ray_dir;               /* Normalized ray direction     */
  573.   double a,b;
  574.                     
  575.   /* Calculate delta form-factor from each sample */
  576.   ff = 0.0;
  577.   for (i=0;i<num_samples;i++) {
  578.     ray_dir = *vsub(&(samples[i].pos), &v);
  579.     a = vlen(vsub(&samples[i].pos, 
  580.           vmiddle(&samples[i].corners[2], &samples[i].corners[1])));
  581.     if (p->class == TRIANGLE)
  582.       b = vlen(vsub(&samples[i].pos, &samples[i].corners[2]));
  583.     else
  584.       b = vlen(vsub(&samples[i].pos, 
  585.             vmiddle(&samples[i].corners[3], &samples[i].corners[2])));
  586.  
  587.     norm_ray_dir = ray_dir;
  588.     norm(&norm_ray_dir);
  589.     w = samples[i].area;
  590.     if (shoot_ray(v, ray_dir, samples[i], &t)) {
  591.       cos_angle_v = dot(&n, &norm_ray_dir);
  592.       cos_angle_p = dot(&(samples[i].norm), vnegate(&norm_ray_dir));
  593.       if (cos_angle_p == 0.0) cos_angle_p = cos(DTOR * 90.0);
  594.  
  595.       if (a == b) /* Is a circle */
  596.     dff = (cos_angle_v * cos_angle_p) / ((M_PI * t*t) + samples[i].area);
  597.       else       /* Is an ellipse */
  598.     dff = (cos_angle_v * cos_angle_p) / 
  599.       (M_PI * sqrt( (t*t+a*a) * (t*t+b*b) ) );
  600.       dff = dff * w;
  601.  
  602.       /* Take area weighted average of samples as form factor */
  603.       ff += dff;
  604.     }
  605.   }
  606.   return(ff);
  607. }
  608.     
  609. /**********************************************************************/
  610. /* Find form factors between all element vertices of receiver and all */
  611. /* source samples.                                                    */
  612. /**********************************************************************/
  613. double ff_patches(src,rec,ff_rs_vtx,
  614.           num_src_samples, src_samples,
  615.           corners)
  616.      Polygon *src;                         /* Source */
  617.      Polygon *rec;                         /* Receiver */
  618.      double ff_rs_vtx[MAX_PATCH_VTX];      /* Receiver form factors */
  619.      int num_src_samples;                  /* Number of source samples */
  620.      Sampletype src_samples[MAX_SAMPLES];  /* Source samples */
  621.      Vector corners[MAX_PATCH_VTX];        /* Corners of source */
  622. {
  623.   int i;
  624.   double ff_rs;                       /* Form factor total */
  625.   double element_ff;
  626.   Vector element_pos, element_norm;
  627.   Vector rec_corners[MAX_PATCH_VTX];    /* Samples for analytic */
  628.   Sampletype rec_samples[MAX_SAMPLES];
  629.   int num_rec_samples;  
  630.  
  631.   float ff_start, ff_end;
  632.   
  633.   ff_rs = 0.0;
  634.   ff_start = Cpu_Time(&tstats.utime, &tstats.stime);  
  635.   
  636.   if (FF_Options.fftype == NUMERICAL_FF) {
  637.     
  638.     /* Check from all corners of receiver */
  639.     for(i=0;i<rec->numVert;i++) { 
  640.       ff_rs_vtx[i] = 0.0;
  641.       element_pos = rec->vert[i]->pos;
  642.       element_norm = rec->normal[0];
  643.       
  644.       /* Check perpendicular condition here. If so sample from 
  645.      centers as well. N/F */
  646.       /* if (dot(&element_norm, &src_samples[0].norm) == 0) { */
  647.  
  648.       if (FF_Options.sample_shape == DISC_FF)
  649.     element_ff = ff_da_disc(element_pos, element_norm,
  650.                 src, src_samples, num_src_samples);
  651.       else
  652.     element_ff = ff_da_ellipse(element_pos, element_norm,
  653.                 src, src_samples, num_src_samples);    
  654.       if (element_ff > 0.0) {
  655.     ff_rs += element_ff;                /* Sum total ff */
  656.     ff_rs_vtx[i] = element_ff;          /* Store vertex ff */
  657.       }
  658.     } 
  659.   } 
  660.  
  661.   else if (FF_Options.fftype == ANALYTIC_FF) {
  662.  
  663.     /* Check from centres of sample areas on receiver */
  664.     for(i=0;i<rec->numVert;i++) 
  665.       rec_corners[i] = rec->vert[i]->pos;
  666.     num_rec_samples = ff_Generate_Patch_Samples(rec,rec_samples, rec_corners, 
  667.                         rec->area, rec->normal[0]);
  668.     for (i=0;i<num_rec_samples;i++) { 
  669.       ff_rs_vtx[i] = ff_da_analytic(src, &rec_samples[i].pos, corners,
  670.                     &rec->normal[0]);
  671.       ff_rs_vtx[i] /= (double)num_rec_samples;
  672.       ff_rs += ff_rs_vtx[i];
  673.     }
  674.   }
  675.  
  676.   ff_end = Cpu_Time(&tstats.utime, &tstats.stime);  
  677.   tstats.avg_ff += (ff_end - ff_start);
  678.  
  679.   if (Option.debug) {
  680.     printf("\t+ Receiver %s%d to source %s%d ff: %g\n\t  with vertex ff: ", 
  681.        rec->name, rec->id, src->name, src->id,ff_rs);
  682.     for(i=0;i<rec->numVert;i++) printf("%g ",ff_rs_vtx[i]);
  683.     printf("\n");
  684.   }
  685.   return (ff_rs);
  686. }
  687.  
  688. /*====================================================================*/
  689. /* Code for form-factor refine from [Hanrahan91]. Currently not used. */
  690. /*====================================================================*/
  691. /**********************************************************************/
  692. /* Create new polygon for poly list of poly, update back of list */
  693. /**********************************************************************/
  694. void LinkPoly_ToPoly(to_poly, link_poly, ff_to_link)
  695.      Polygon *to_poly, *link_poly;
  696.      double ff_to_link;
  697. {
  698.   PolyList *pl;
  699.  
  700.   pl = (PolyList *) malloc(sizeof(PolyList));
  701.   pl->patch = link_poly;
  702.   pl->next = (PolyList *)0;
  703.  
  704.   if (to_poly->polyhead.num_polys == 0) { /* No links with any polys yet */
  705.     to_poly->Links = pl;
  706.     to_poly->polyhead.front = to_poly->polyhead.back = to_poly->Links;
  707.   } else {          /* Add to end of list */
  708.     to_poly->polyhead.back->next = pl;
  709.     to_poly->polyhead.back = pl;
  710.   }
  711.   to_poly->polyhead.num_polys++;
  712. }
  713.  
  714. /**********************************************************************/
  715. /* Polygon p can "see" q and q can "see" p                            */
  716. /**********************************************************************/
  717. void ff_Link_Patches(p,q, ff_pq, ff_qp)
  718.      Polygon *p, *q;
  719.      double ff_pq, ff_qp;
  720. {
  721.   if (Option.debug)
  722.     printf("\tLinking polygon %s%d and poly %s%d.\n",
  723.        p->name, p->id, q->name, q->id);
  724.   LinkPoly_ToPoly(p,q, ff_pq);
  725.   LinkPoly_ToPoly(q,p, ff_qp);
  726. }
  727.  
  728. /**********************************************************************/
  729. /* Perform BF refinement on 2 patches                                 */
  730. /**********************************************************************/
  731. void ff_BF_Refine(src,rec)
  732.      Polygon *src, *rec;
  733. {
  734.   int i;
  735.   double ff_rs, ff_sr;
  736.   double ff_difference;
  737.   double subdiv_src_area, subdiv_rec_area;
  738.   double ff_rs_vtx[MAX_PATCH_VTX];
  739.   double ff_sr_vtx[MAX_PATCH_VTX];
  740.  
  741.   int num_src_samples;          
  742.   int num_rec_samples;          
  743.   Sampletype src_samples[MAX_SAMPLES];
  744.   Sampletype rec_samples[MAX_SAMPLES];
  745.   Vector src_corners[MAX_PATCH_VTX];
  746.   Vector rec_corners[MAX_PATCH_VTX];
  747.   
  748.   /* Find corners of source and receiver, and generate sample 
  749.      points on source and receiver patch */
  750.   for(i=0;i<src->numVert;i++) src_corners[i] = src->vert[i]->pos;
  751.   num_src_samples = ff_Generate_Patch_Samples(src,src_samples, 
  752.                           src_corners, src->area, 
  753.                           src->normal[0]);
  754.   for(i=0;i<rec->numVert;i++) rec_corners[i] = rec->vert[i]->pos;
  755.   num_rec_samples = ff_Generate_Patch_Samples(rec,rec_samples, 
  756.                           rec_corners, rec->area, 
  757.                           rec->normal[0]);
  758.  
  759.   /* Find form factors from source to receiver and vice versa */
  760.   ff_sr = ff_patches(src,rec,ff_sr_vtx, num_src_samples,
  761.              src_samples, src_corners);
  762.   ff_rs = ff_patches(rec,src,ff_rs_vtx, num_rec_samples,
  763.              rec_samples, rec_corners);
  764.   ff_difference = fabs(ff_sr - ff_rs);
  765.   subdiv_src_area = src->area / 4.0;
  766.   subdiv_rec_area = rec->area / 4.0;
  767.  
  768.   if ((ff_rs < FF_Options.F_min) && (ff_sr < FF_Options.F_min))
  769.     ff_Link_Patches(src,rec, ff_sr, ff_rs);
  770.   else if (ff_difference > FF_Options.F_diff_poly) {
  771.  
  772.     if (ff_sr > ff_rs) {
  773.       if (subdiv_rec_area > FF_Options.min_element_area) {
  774.     /* Subdivide polygon rec and refine */
  775.     Subdivide_Polygon(rec);
  776.     for (i=0;i<MAX_PATCH_CHILDREN;i++) 
  777.       ff_BF_Refine(src,rec->child[i]);
  778.       } else 
  779.     ff_Link_Patches(src,rec, ff_sr, ff_rs);
  780.     } else {
  781.       if (subdiv_src_area > FF_Options.min_element_area) {
  782.     /* Subdivide polygon src and refine */
  783.     Subdivide_Polygon(src);
  784.     for (i=0;i<MAX_PATCH_CHILDREN;i++) 
  785.       ff_BF_Refine(rec,src->child[i]);
  786.       } else 
  787.     ff_Link_Patches(src,rec, ff_sr, ff_rs);
  788.     }
  789.   }
  790. }
  791.  
  792.